home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / blas / zher2.f < prev    next >
Text File  |  1997-06-25  |  9KB  |  253 lines
  1. *
  2. ************************************************************************
  3. *
  4.       SUBROUTINE ZHER2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA )
  5. *     .. Scalar Arguments ..
  6.       COMPLEX*16         ALPHA
  7.       INTEGER            INCX, INCY, LDA, N
  8.       CHARACTER*1        UPLO
  9. *     .. Array Arguments ..
  10.       COMPLEX*16         A( LDA, * ), X( * ), Y( * )
  11. *     ..
  12. *
  13. *  Purpose
  14. *  =======
  15. *
  16. *  ZHER2  performs the hermitian rank 2 operation
  17. *
  18. *     A := alpha*x*conjg( y' ) + conjg( alpha )*y*conjg( x' ) + A,
  19. *
  20. *  where alpha is a scalar, x and y are n element vectors and A is an n
  21. *  by n hermitian matrix.
  22. *
  23. *  Parameters
  24. *  ==========
  25. *
  26. *  UPLO   - CHARACTER*1.
  27. *           On entry, UPLO specifies whether the upper or lower
  28. *           triangular part of the array A is to be referenced as
  29. *           follows:
  30. *
  31. *              UPLO = 'U' or 'u'   Only the upper triangular part of A
  32. *                                  is to be referenced.
  33. *
  34. *              UPLO = 'L' or 'l'   Only the lower triangular part of A
  35. *                                  is to be referenced.
  36. *
  37. *           Unchanged on exit.
  38. *
  39. *  N      - INTEGER.
  40. *           On entry, N specifies the order of the matrix A.
  41. *           N must be at least zero.
  42. *           Unchanged on exit.
  43. *
  44. *  ALPHA  - COMPLEX*16      .
  45. *           On entry, ALPHA specifies the scalar alpha.
  46. *           Unchanged on exit.
  47. *
  48. *  X      - COMPLEX*16       array of dimension at least
  49. *           ( 1 + ( n - 1 )*abs( INCX ) ).
  50. *           Before entry, the incremented array X must contain the n
  51. *           element vector x.
  52. *           Unchanged on exit.
  53. *
  54. *  INCX   - INTEGER.
  55. *           On entry, INCX specifies the increment for the elements of
  56. *           X. INCX must not be zero.
  57. *           Unchanged on exit.
  58. *
  59. *  Y      - COMPLEX*16       array of dimension at least
  60. *           ( 1 + ( n - 1 )*abs( INCY ) ).
  61. *           Before entry, the incremented array Y must contain the n
  62. *           element vector y.
  63. *           Unchanged on exit.
  64. *
  65. *  INCY   - INTEGER.
  66. *           On entry, INCY specifies the increment for the elements of
  67. *           Y. INCY must not be zero.
  68. *           Unchanged on exit.
  69. *
  70. *  A      - COMPLEX*16       array of DIMENSION ( LDA, n ).
  71. *           Before entry with  UPLO = 'U' or 'u', the leading n by n
  72. *           upper triangular part of the array A must contain the upper
  73. *           triangular part of the hermitian matrix and the strictly
  74. *           lower triangular part of A is not referenced. On exit, the
  75. *           upper triangular part of the array A is overwritten by the
  76. *           upper triangular part of the updated matrix.
  77. *           Before entry with UPLO = 'L' or 'l', the leading n by n
  78. *           lower triangular part of the array A must contain the lower
  79. *           triangular part of the hermitian matrix and the strictly
  80. *           upper triangular part of A is not referenced. On exit, the
  81. *           lower triangular part of the array A is overwritten by the
  82. *           lower triangular part of the updated matrix.
  83. *           Note that the imaginary parts of the diagonal elements need
  84. *           not be set, they are assumed to be zero, and on exit they
  85. *           are set to zero.
  86. *
  87. *  LDA    - INTEGER.
  88. *           On entry, LDA specifies the first dimension of A as declared
  89. *           in the calling (sub) program. LDA must be at least
  90. *           max( 1, n ).
  91. *           Unchanged on exit.
  92. *
  93. *
  94. *  Level 2 Blas routine.
  95. *
  96. *  -- Written on 22-October-1986.
  97. *     Jack Dongarra, Argonne National Lab.
  98. *     Jeremy Du Croz, Nag Central Office.
  99. *     Sven Hammarling, Nag Central Office.
  100. *     Richard Hanson, Sandia National Labs.
  101. *
  102. *
  103. *     .. Parameters ..
  104.       COMPLEX*16         ZERO
  105.       PARAMETER        ( ZERO = ( 0.0D+0, 0.0D+0 ) )
  106. *     .. Local Scalars ..
  107.       COMPLEX*16         TEMP1, TEMP2
  108.       INTEGER            I, INFO, IX, IY, J, JX, JY, KX, KY
  109. *     .. External Functions ..
  110.       LOGICAL            LSAME
  111.       EXTERNAL           LSAME
  112. *     .. External Subroutines ..
  113.       EXTERNAL           XERBLA
  114. *     .. Intrinsic Functions ..
  115.       INTRINSIC          DCONJG, MAX, DBLE
  116. *     ..
  117. *     .. Executable Statements ..
  118. *
  119. *     Test the input parameters.
  120. *
  121.       INFO = 0
  122.       IF     ( .NOT.LSAME( UPLO, 'U' ).AND.
  123.      $         .NOT.LSAME( UPLO, 'L' )      )THEN
  124.          INFO = 1
  125.       ELSE IF( N.LT.0 )THEN
  126.          INFO = 2
  127.       ELSE IF( INCX.EQ.0 )THEN
  128.          INFO = 5
  129.       ELSE IF( INCY.EQ.0 )THEN
  130.          INFO = 7
  131.       ELSE IF( LDA.LT.MAX( 1, N ) )THEN
  132.          INFO = 9
  133.       END IF
  134.       IF( INFO.NE.0 )THEN
  135.          CALL XERBLA( 'ZHER2 ', INFO )
  136.          RETURN
  137.       END IF
  138. *
  139. *     Quick return if possible.
  140. *
  141.       IF( ( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) )
  142.      $   RETURN
  143. *
  144. *     Set up the start points in X and Y if the increments are not both
  145. *     unity.
  146. *
  147.       IF( ( INCX.NE.1 ).OR.( INCY.NE.1 ) )THEN
  148.          IF( INCX.GT.0 )THEN
  149.             KX = 1
  150.          ELSE
  151.             KX = 1 - ( N - 1 )*INCX
  152.          END IF
  153.          IF( INCY.GT.0 )THEN
  154.             KY = 1
  155.          ELSE
  156.             KY = 1 - ( N - 1 )*INCY
  157.          END IF
  158.          JX = KX
  159.          JY = KY
  160.       END IF
  161. *
  162. *     Start the operations. In this version the elements of A are
  163. *     accessed sequentially with one pass through the triangular part
  164. *     of A.
  165. *
  166.       IF( LSAME( UPLO, 'U' ) )THEN
  167. *
  168. *        Form  A  when A is stored in the upper triangle.
  169. *
  170.          IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
  171.             DO 20, J = 1, N
  172.                IF( ( X( J ).NE.ZERO ).OR.( Y( J ).NE.ZERO ) )THEN
  173.                   TEMP1 = ALPHA*DCONJG( Y( J ) )
  174.                   TEMP2 = DCONJG( ALPHA*X( J ) )
  175.                   DO 10, I = 1, J - 1
  176.                      A( I, J ) = A( I, J ) + X( I )*TEMP1 + Y( I )*TEMP2
  177.    10             CONTINUE
  178.                   A( J, J ) = DBLE( A( J, J ) ) +
  179.      $                        DBLE( X( J )*TEMP1 + Y( J )*TEMP2 )
  180.                ELSE
  181.                   A( J, J ) = DBLE( A( J, J ) )
  182.                END IF
  183.    20       CONTINUE
  184.          ELSE
  185.             DO 40, J = 1, N
  186.                IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN
  187.                   TEMP1 = ALPHA*DCONJG( Y( JY ) )
  188.                   TEMP2 = DCONJG( ALPHA*X( JX ) )
  189.                   IX    = KX
  190.                   IY    = KY
  191.                   DO 30, I = 1, J - 1
  192.                      A( I, J ) = A( I, J ) + X( IX )*TEMP1
  193.      $                                     + Y( IY )*TEMP2
  194.                      IX        = IX        + INCX
  195.                      IY        = IY        + INCY
  196.    30             CONTINUE
  197.                   A( J, J ) = DBLE( A( J, J ) ) +
  198.      $                        DBLE( X( JX )*TEMP1 + Y( JY )*TEMP2 )
  199.                ELSE
  200.                   A( J, J ) = DBLE( A( J, J ) )
  201.                END IF
  202.                JX = JX + INCX
  203.                JY = JY + INCY
  204.    40       CONTINUE
  205.          END IF
  206.       ELSE
  207. *
  208. *        Form  A  when A is stored in the lower triangle.
  209. *
  210.          IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
  211.             DO 60, J = 1, N
  212.                IF( ( X( J ).NE.ZERO ).OR.( Y( J ).NE.ZERO ) )THEN
  213.                   TEMP1     = ALPHA*DCONJG( Y( J ) )
  214.                   TEMP2     = DCONJG( ALPHA*X( J ) )
  215.                   A( J, J ) = DBLE( A( J, J ) ) +
  216.      $                        DBLE( X( J )*TEMP1 + Y( J )*TEMP2 )
  217.                   DO 50, I = J + 1, N
  218.                      A( I, J ) = A( I, J ) + X( I )*TEMP1 + Y( I )*TEMP2
  219.    50             CONTINUE
  220.                ELSE
  221.                   A( J, J ) = DBLE( A( J, J ) )
  222.                END IF
  223.    60       CONTINUE
  224.          ELSE
  225.             DO 80, J = 1, N
  226.                IF( ( X( JX ).NE.ZERO ).OR.( Y( JY ).NE.ZERO ) )THEN
  227.                   TEMP1     = ALPHA*DCONJG( Y( JY ) )
  228.                   TEMP2     = DCONJG( ALPHA*X( JX ) )
  229.                   A( J, J ) = DBLE( A( J, J ) ) +
  230.      $                        DBLE( X( JX )*TEMP1 + Y( JY )*TEMP2 )
  231.                   IX        = JX
  232.                   IY        = JY
  233.                   DO 70, I = J + 1, N
  234.                      IX        = IX        + INCX
  235.                      IY        = IY        + INCY
  236.                      A( I, J ) = A( I, J ) + X( IX )*TEMP1
  237.      $                                     + Y( IY )*TEMP2
  238.    70             CONTINUE
  239.                ELSE
  240.                   A( J, J ) = DBLE( A( J, J ) )
  241.                END IF
  242.                JX = JX + INCX
  243.                JY = JY + INCY
  244.    80       CONTINUE
  245.          END IF
  246.       END IF
  247. *
  248.       RETURN
  249. *
  250. *     End of ZHER2 .
  251. *
  252.       END
  253.